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Abstract: We show how to efficiently calculate the signal in optical 
coherence tomography (OCT) systems due to the ballistic photons, the 
quasi-ballistic photons, and the photons that undergo multiple diffusive 
scattering using Monte Carlo simulations with importance sampling. This 
method enables the calculation of these three components of the OCT signal 
with less than one hundredth of the computational time required by the 
conventional Monte Carlo method. Therefore, it can be used as a design tool 
to characterize the performance of OCT systems, and can also be used in the 
development of novel signal processing techniques that can extend the 
imaging range of OCT systems. We investigate the parameter dependence 
of our importance sampling method and we validate it by comparison to an 
existing method. 
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1. Introduction 

Optical coherence tomography (OCT) is a sub-surface imaging technique that produces tissue 
images one to two orders of magnitude higher resolution than ultrasound imaging and an 
imaging depth that can reach up to 3 mm, depending on the optical properties of the tissue [1], 
and can be integrated with optical fiber probes and also allow spectroscopic characterization 
of tissue [2]. Computer simulations of light transport in multi-layered turbid media are an 
effective way to theoretically investigate light transport in organic tissues [1], which can also 
be applied in the development of optical coherence tomography (OCT) systems [2]. 

Yao and Wang [3] proposed an importance sampling implementation to speed up Monte 
Carlo simulations in time-domain OCT (TD-OCT) systems. In [4], we reviewed previous 
methods proposed to speed up Monte Carlo simulations and we showed how one can extend 
the method in [3] using importance sampling to significantly speed up the calculation of the 
OCT signal from the entire depth range imaged that results from the ballistic and the quasi- 
ballistic scattered photons, the Class I reflectance, without significant residual bias in the 
statistical result. That importance sampling method consists of applying multiple biases 
towards the source after the first biased scattering is applied to the photon packets, and using a 
photon splitting procedure. The Class I reflectance is produced by the photons that only 
undergo ballistic or quasi-ballistic propagation, besides undergoing back-scattering at the 
layer being imaged. Importance sampling, tailored to each particular application, is a 
statistical technique that has also been applied in optical communications [5-7], confocal 
microscopy [8], atmospheric optics [9], and in diffuse optical tomography [10]. 

In this paper we show that the implementation used in [4] can be enhanced to also 
efficiently calculate the OCT signal due to the detected photons that undergo multiple 
diffusive scattering in tissue, the Class II reflectance. The Class II reflectance interferes with 
the Class I reflectance and, as a consequence, limits the ability of OCT systems to image 
deeper into the tissue [11]. The importance sampling method in [4] consists of randomly 
applying multiple biases from the first biased scattering towards the tip of the collecting 
optical fiber. In the method proposed here, the application of additional biased scatterings 
towards the collecting optics, after the first biased scattering, are not carried out automatically 
as in the previous method. Instead, the bias is applied with a probability p specified for the 
ensemble of Monte Carlo simulations [12]. In the event in which a bias towards the collecting 
optics is not applied, the scattering will follow the unbiased model. Moreover, the biased 
scatterings in this method can only produce scatterings in the half-sphere along the apparent 
direction of the collecting optics, as opposed to scattering in all directions as in the previous 
method. These procedures enable a significant increase in the number of collected photon that 
undergo multiple diffusive scattering, along with a decrease in the range of values of the 
likelihood ratio that accelerate the convergence of the calculation of the Class II reflectance. 
The combination of these two factors leads to faster convergence of the Monte Carlo 
simulations for the calculation of both the Class I and the Class II reflectances. 

In Sec. 2, we described the variance reduction method that we developed to speed up the 
calculation of the reflectance using Monte Carlo simulations. In Sec. 3, we show numerical 
results for the standard Monte Carlo method and for our importance sampling 
implementations with different parameters. In that section, we validated the importance 
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sampling method by comparing it against extensive standard Monte Carlo simulations, and we 
demonstrate the effectiveness of our importance sampling implementation. 

2. Importance sampling applied to OCT 

The reflectance in turbid media, such as organic tissue, obtained with OCT systems can be 
easily calculated by modeling light as electromagnetic fields using the first Born 
approximation [13], without which computer simulations using the full field model would not 
be practical. However, this approximation only includes the contribution of the ballistic 
photons. Therefore, this simplified model cannot account for the photons that undergo 
multiple diffusive scattering in tissue, which can limit the imaging depth of OCT systems. 
One alternative to this method is the use of the Monte Carlo method do model light transport 
in tissue, in which the light is modeled as photon particles that do not interact with one 
another. This method is applicable when the coherence length of the light is short when 
compared to both the mean path length in tissue and the widths of the tissue layers, which 
corresponds to most practical cases in which OCT imaging is exploited. 

We implemented our importance sampling method for Monte Carlo modeling of light 
transport in multi-layered tissues (MCML) by modifying a C-language software package that 
is available for download from the web site of the Oregon Medical Laser Center [14,15]. 
MCML can be used to simulate an ensemble of photon packets launched in a steady-state one- 
dimensional beam, normal to the surface of the topmost tissue layer. The scattering events that 
arise from the simulated tissue are characterized by two random angles that determine the 
future direction of the photon packet in three-dimensional space after the scattering. As in the 
MCML, to account for the photon packet scattering with arbitrary anisotropy factor g of the 
tissue, in the unbiased scatterings we use the Henyey-Greenstein probability density function 
that is defined as 



where 6 S is the angle between the photon packet propagation direction u prior to the 
scattering and the new scattering direction u ' . After rotating away from the previous 
propagation direction u by the angle 0 s , so that cos 0 s = u • u ' , the scattering direction u ? is 

rotated around u by an angle (j) that is randomly picked from a uniform probability density 
function from 0 to 2n. 

2.1. First biased scattering 

Considering that the probability that a photon packet is scattered towards the apparent position 
of the collecting optics at any given layer is very low, since the anisotropy factor g of tissue is 
close to 1 , and the probability of reflection also decreases with the depth, the convergence of 
the calculation of the Class I and the Class II reflectances using standard Monte Carlo 
simulations is very slow, as shown in [3]. In [4] we proposed an importance sampling 
implementation for the Class I reflectance that was based on biased scattering towards the 
apparent position of the collecting optics and a photon splitting procedure followed by 
successive biased scatterings towards the apparent position of the collecting optics, whose 
direction we defined as the unit vector v . Because that method was not as efficient in the 
calculation of the Class II reflectance, we significantly enhanced that method by 
implementing two modifications. The first modification consists of using a biased probability 
density function for the first biased scattering that produces random scattering with an angle 
no lesser than 90 degrees away from the direction to the apparent position of the collecting 
optics. This biased probability density function, which was based on the Henyey-Greenstein 
probability density function in Eq. (1), is given by 
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(1) 



f B (cos0 B ) = 
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0, otherwise 



where a is a bias coefficient in the range (0,1). After randomly picking a biased angle 0 B 
away from the direction of the apparent position of the collecting optics, the biased 
direction v , so that cos 9 B = v • u ' , the resultant biased scattering direction u ' is rotated 

around v by an angle $ that is randomly picked from a uniform probability density function 
from 0 to 27i. This last procedure is equivalent to the one used in the MCML software package 
to enable a full three-dimensional scattering. Then, the scattered photon packet is associated 
with a quantity that is defined as the likelihood ratio [5-7], which ensure converge of the 
statistical result towards the expected value. The likelihood ratio of the photon packet using 
the biased probability density function in Eq. (2) is given by 



/ x fur (COS I- 
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5 J 



where cos# 5 =u-u' is a function of cos# g , which is statistically drawn from the probability 

density function in Eq. (2) that is used to define the new propagation direction u' of the 
photon packet at the first biased scattering. A schematic representation of the angles and 
vectors used in the biased and in the unbiased scatterings is shown in Fig. 1 of [4]. 

Other probability density functions can also effectively speed up the calculation using this 
method, provided that they significantly increase the probability that the photon packet is 
scattered towards the apparent position of the collecting optics. 

2.2. Additional biased scatterings 

The second enhancement that we made to the variance reduction method in [4] consists of 
applying additional biased scatterings towards the apparent direction of the collecting fiber v 
with probability 0 < p < 1, as opposed to p = 1 as it is the previous method. If a biased 
scattering is not applied at a given point in the tissue, the photon packet will undergo an 
unbiased scattering at that location according to the probability density function in Eq. (1). 
The likelihood ratio of this scattering, whether the biased or the unbiased probability density 
function is randomly selected, is calculated according to the expression 

L(cos0 B ) = - t ^ s - j — i r. (4) 

V B) pf B (cos0 B ) + (l-p)f HG (cos0 s ) 

If the biased function f B (cos# g ) is selected to draw a random value of cos# g , which is 

an event with probability p, cos 9 s = u • u ' is a function of cos 6 B that is statistically drawn 

from the probability density function in Eq. (2). Otherwise, which is a complementary event 

with probability 1 - p, the unbiased function f HG (cos# 5 ) is selected to draw a random value 

of cos# s and cos# g =v-u' is a function of cos# s . Because each random scattering is 

independent from the other scatterings, the likelihood ratio of each photon packet is equal to 
the product of all the likelihood ratios of all the biased scatterings of that photon packet. 

2.3. Calculation of the reflectance 

The Class I and the Class II reflectances at depth z are obtained by calculating the mean value 
of the indicator functions I x and I 2 of the spatial and temporal filter of the Class I reflectance 
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and the Class II reflectance, respectively, for all the photon packets (samples) in the ensemble. 
The indicator function I x and I 2 of the spatial and temporal filter for the zth photon packet is 
defined as 

[l, / >|As. -2z \,r.<d ,0.<0 , |As. -2z|< / 

r / -\ ' c \ i max ' i max' z,i max' / c , c \ 

7 i( z >0H ' . ( 5 ) 



0, otherwise 



and 



fl, / <|As. -2z \,r.<d ,0 <6 , lAs*. -2z\<1 

r / -\ ' c max ' i max' z,i max' c / s\ 

I 2 (z,i) = \ 1 m (6) 

[0, otherwise 

where z is the depth imaged, l c is the coherence length of the source, r t is the distance of the 
zth reflected photon packet to the origin along the plane z = 0, where the collecting optical 
system is located, d max and 0 max are the maximum collecting diameter and angle, respectively, 
6 zi is the angle with the z-axis, which is the axis normal to the tissue interface, As t is the 
optical path, z max is the maximum depth reached by the photon packet. A detected photon 
packet is considered a Class II photon packet if it does not reach a depth that is consistent with 
its optical path, so that it interferes constructively with corresponding detected Class I photons 
packets without bringing any information from the depth in which those Class I photons 
packets were reflected. For simplicity, the indicator functions in Eq. (5) and in Eq. (6) were 
defined with a square time gating as in [3]. The estimated values of the Class I and Class II 
reflectances and their respective variances are given by the following expressions 

R i,i ( z ) = 77 Z A ( z > i)L(iW(j) (7) 

and 

o£ 2 (z) = N (x-l) tt^ (Z ' i)L(iW(i) ~ Rl2 (Z) ] ' (8) 

where W(i) is the weight of the z'th photon packet in MCML, which is a quantity affected by 
the absorption coefficient at the scattering points, and L(7) is the product of the likelihood 
rations of all the biased scatterings that affected the z'th photon packet. Using the Monte Carlo 
method with the importance sampling implementation described in this section, the 
calculation of the reflectances in Eq. (7) converge two to three orders of magnitude faster with 
the number of samples TV than the standard Monte Carlo method used in MCML. 

2.4. Generation of random biased angles 

To generate random angles according to the biased probability density function in Eq. (2) we 
use the pseudo-random number generator of the Gnu Scientific Library [16]. That random 
number generator produces pseudo-random numbers uniformly distributed from 0 to 1 , which 
we convert to the probability density function in Eq. (2) according to the following conversion 
formula 



1 | , 
cos@ Bi = — \a +1- 
2a 



1 1 



Vtf 2 +1 ) Va 2 +1 



1 



(9) 



where u t is sampled from the random number generator of the Gnu Scientific Library. 
Equation (9) was derived based on the probability theory in [17]. 
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3. Numerical results 



We validate the importance sampling method for Monte Carlo simulations that we developed 
by comparing it against extensive standard Monte Carlo simulations. We simulate a turbid 
media that consists of multiple layers, shown schematically in Fig. 1 . The tissue extends from 
0 to 1 mm, consisting primarily of a turbid layer with absorption coefficient ju a = 1.5 cm" 1 and 
a scattering coefficient ju s = 60 cm" 1 , but also contains five thin layers with absorption 
coefficient ju a = 3 cm" 1 and a scattering coefficient ju s = 120 cm" 1 . These five thin layers with 
higher scattering coefficient are located from 200 jum to 215 jum, from 365 jum to 395 jum, 
from 645 jum to 660 jum, from 760 jum to 775 jum, and from 900 jum to 915 jum. We assume 
that this tissue has the same refractive index n = 1 and an anisotropy factor g = 0.9, as in [3]. 
In [4] we showed that our method is robust in the presence of refractive index mismatch along 
the tissue, in which the apparent position of the collecting optics is different in each layer due 
to refraction at the interfaces. We simulate a TD-OCT system that is delivered and collected 
by the tip of an optical fiber with a radius of 10 jum and an acceptance angle of 5 degrees. For 
simplicity, the light source is assumed to be a one-dimensional light beam propagating along 
the vertical direction as in [3,4,10], since the purpose of this study is to validate and 
demonstrate the effectiveness of the importance sampling implementation that we developed. 



OPTICAL 
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y-axis 
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X-axis 
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Fig. 1. Schematic representation of a simulation setup similar to [5]. 

In Figs. 2 and 3, we show results obtained with 10 8 Monte Carlo simulations with 
importance sampling, which has a computational cost of about 9 x 1 0 8 standard Monte Carlo 
simulations in this particular simulation setup because of the computational cost of the photon 
splitting procedure [4]. The computational cost increase of this importance sampling method 
depends on the target depth range and on the photon mean free path in the tissue. The target 
depth range of these simulations was set from 0 mm to 1 mm. Therefore, every single photon 
scattering that occurs in the depth range from 0 mm to 1 mm is biased. We run the Monte 
Carlo simulations with importance sampling with the bias coefficient a = 0.925 and additional 
bias probability p = 0.5. The results shown in Figs. 2 and 3 indicate that our new importance 
sampling procedure reduces the computational cost to obtain the Class I diffuse reflectance by 
about three orders of magnitude when compared to the standard Monte Carlo method. 

We used a computer with the AMD Opteron 246 processor with a clock of 2 GHz and 
4GB of RAM to run all the simulations presented. The simulation using the standard Monte 
Carlo method, whose results are shown in Figs. 2 and 3, required eight days, 22 hours, and 7 
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minutes of computer time, while the simulation using the importance sampling method for 
Monte Carlo required only 1 hour and 53 minutes of computer time. We observed that the 
confidence intervals of the results obtained using the standard Monte Carlo method are 
significantly larger than the ones obtained with importance sampling for the results shown in 
Fig. 3, even though the standard Monte Carlo simulations required 113 times the 
computational time of the simulations with importance sampling. It would have been 
necessary to increase the number of samples in the standard Monte Carlo simulations by one 
order of magnitude (89 days of computer simulation) to obtain confidence intervals of the 
Class I reflectance comparable to those obtained using importance sampling. In Fig. 3, we 
also observed that this method reduces the computational cost of calculating the Class II 
reflectance by more than two orders of magnitude. 

In Fig. 4 we show the dependence of the relative error of the calculation of the Class I and 
the Class II reflectances at two different depths: 400 jum and 670 jum with respect to the bias 
coefficient a for/? = 0.5. The depths at 400 jum and 670 urn correspond to the tissue regions 




Depth (mm) 



Fig. 2. The Class I reflectance, shown with thick solid black curve, and the Class II reflectance, 
shown with thin solid red curve, as a function of the depth for the importance sampling 
implementation described in Sec. 2 with 10 8 samples. The pink short dashed curve and the blue 
long dashed curve are results of 10 11 standard Monte Carlo simulations of the Class I 
reflectance and the Class II reflectance, respectively. 



10" 




0.64 



Depth (mm) 



0.68 



Fig. 3. The reflectance results shown in Fig. 2 for the depth interval from 640 urn to 680 um. 
The error bars shown for every other point were estimated in the same ensemble of simulations. 
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Fig. 4. The relative error in the calculation of the reflectance using importance sampling as a 
function of the bias coefficient a at 400 urn and at 670 urn of depth for p = 0.5. 




Q I I I I I I I I I L 

^ Probability of Bias 



Fig. 5. The relative error in the calculation of the reflectance using importance sampling as a 
function of the probability of additional bias p at 400 um and at 670 um of depth for a = 0.9. 

near the second and the third peaks of the reflectance. The relative error is defined as the ratio 
between the standard deviation, which is the square root of the variance in Eq. (8), divided by 
the estimated value of the reflectance in Eq. (7). The Class I reflectance has its minimum 
relative error at 400 jum close to a = 0.925, but the minimum shifts to about a = 0.95 um at 
670 jum because a stronger bias is necessary to collect more samples from deeper regions in 
the tissue. However, as the bias coefficient is increased towards 1, larger variations in the 
likelihood ratio due to very strong bias leads to an increase in the relative error with the bias 
coefficient. The Class II reflectance has its minimum relative error at 400 um close to a = 
0.91, and shifts to only about a = 0.925 jum at 670 um. That is lower than the optimum bias 
coefficient observed in the Class I reflectance because strong bias leads to an increase in the 
number of ballistic and quasi-ballistic photons and a decrease in the number of collected 
photons that were scattered multiple times in the tissue. Figure 4 shows that there is a region 
between 0.9 and 0.95 for the bias parameter a that enables a rapid convergence of the 
calculation of the Class I and Class II reflectance with the Monte Carlo method with this 
importance sampling implementation. 

In Fig. 5 we show the dependence of the relative error of the calculation of the Class I and 
the Class II reflectances at two different depths: 400 um and 670 um, with respect to the 
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probability of additional bias p for a = 0.9. At the depth of 400 jum, the relative error is 
minimized when p = 0.5 for both the Class I and the Class II reflectances. At the depth of 670 
jum the relative error of the Class I reflectance is minimized at p = 0.6, while the Class II 
reflectance has its relative error minimized at p = 0.55. Nevertheless, the relative error of the 
numerical results is not very sensitive to the value of p when p is between 0.3 and 0.7. At p = 
1 the bias used is comparable to the one presented in [4], in which all the additional 
scatterings (after the first biased scattering) are biased. In that case, the Class II reflectance at 
670 jum of depth has a relative error that is nearly twice the minimum value obtained when p 
= 0.55, which implies into a convergence in the new method that is four times faster than the 
previous method. At p = 0, which corresponds to the parameter regime in which only one bias 
scattering is applied, this method is the least effective in the reflectance calculation. 

4. Conclusion 

We developed a new importance sampling method that enables a reduction of the computation 
time of the Class I and the Class II reflectances in TD-OCT by as much as three orders of 
magnitude, which we validated by comparing it against a large ensemble of standard Monte 
Carlo simulations based on MCML. This enables a computation time reduction from several 
days down to tens of minutes. This fast Monte Carlo calculation of TD-OCT signals can be a 
valuable tool in the investigation of the physical process governing both the Class I and the 
Class II reflectances that may enable the development of novel signal processing techniques to 
extend the imaging depth of these systems. Since the photons that undergo multiple 
scatterings in tissue are also responsible for the Class II reflectance in Frequency Domain 
OCT systems, we believe that this importance sampling algorithm can be easily extended to 
those systems. 
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